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1.0  STATEMENT  OF  THE  PROBLEM  STUDIED 


Application  of  holographic  tomography  to  the  flow  around  the  tip 
of  a  revolving  airfoil  was  investigated.  Data  in  the  form  of  interfer- 
ograms  obtained  at  different  angles  obtained  by  AVRADCOM  at  NASA  Ames 
were  used  to  reconstruct  the  three-dimensional  density  field.  The  data 
were  obtained  for  the  flow  around  a  rotating  NACA0012  airfoil  at  a  tip 
mach  number  of  0.9  (See  Figure  1).  Each  interferogram  represented  the 
integrated  map  of  the  optical  path  along  a  given  orientation  relative  to 
the  wing.  A  total  of  forty  interferograms  were  used.  A  typical  inter¬ 
ferogram  is  shown  in  Figure  2.  The  final  objective  of  the  work  was  to 
reconstruct  the  density  field  around  the  tip  of  the  rotating  blade. 

The  work  consisted  of  three  stages: 

a)  Digitization  of  the  interferograms  in  the  phase  maps. 

b)  Reconstruction  of  the  set  of  integrated  phase  maps  into  a 
three-dimensional  map  of  index  of  refraction. 

c)  Development  of  the  density  or  mach  number  maps. 

Each  stage  of  the  work  is  briefly  described  here: 
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Interferogram  Digitization: 

In  the  course  of  this  project,  a  serai-manual  and  an  automated 
technique  for  digitization  of  interferograras  were  developed. 
The  hardware  for  the  manual  digitization  consisted  of  a 
Tektronix  4112  display  terminal  and  a  graphic  tablet  con¬ 
nected  to  a  VAX1 1/780.  An  efficient,  menu-driven  program  for 
digitization  of  the  interferograras  was  constructed.  This 
technique  was  initially  used  for  the  digitization  of  the 
Interferograras . 

An  automated  image  analyzing  system  was  later  developed.  The 
hardware  consisted  of  an  IBM/AT  personal  coaputer,  a  video 
camera  and  a  frame  grabber.  It  had  a  512  x  512  pixel  reso¬ 
lution  which  could  digitize  at  30  frames  per  second  rate. 

The  fringe  digitization  hardware  Is  shown  In  Figure  2. 

The  software  consisted  of  three  n*  i  >x  sets  of  programs.  The 
first  set  of  the  processing  software  Included  most  of  the 
standard  "textbook"  image  processing  procedures'  such  as 
Image  averaging,  filtering,  enhancement,  zooming,  mirror 
imaging  for  the  axlsymmetric  problem  along  with  Sobe 1  or 
Robert  edge  operators.  This  set  of  programs  was  supported  by 
a  commercially  available  Image  processing  software  product. 


The  second  program  set  was  developed  to  determine  fringe 


center-line  coordinates  according  to  the  local  gray  scale 
d istribution.  This  program  was  able  to  detect  fringes  with 


relatively  poor  visibility  and  with  variable  background  gray 
scale. 


The  last  program  was  the  fringe  detection  algorithm  which 
transformed  the  individual  fringe  line  Information  into  a 


R 


phase  map  for  the  entire  pixels.  This  was  achieved  by  sur¬ 
face  spline  Interpolation  technique^,  an  efficient  technique 
to  Interpolate  between  sparse  or  nonuniformly  spaced  data 
points.  It  solves  the  equation  for  elastic  deformation  of  a 
plate  and  assumes  continuity  of  the  function  and  its 
derivatives. 


h)  Tomographic  Reconst  rue  t ion 


A  new  ART  code  based  on  the  iterative  refinement  method  of 
least  squares  solution  for  tomographic  interferometry  was 
developed.  Three  dimensional  field  of  the  refractive  index 
was  reconstructed  from  the  optical  path  data  recorded  on  the 
inter  f erograms . 


The  integral  optical  length  for  negligible  rav  curvature  Is 


expressed  by  the  Radon  transformation^: 
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(1) 


+«• 

l(P,9)  “I  f  n(x,y)  fi[p-xcos(  9)-ysin(  6)  jdxdy 

where  5  is  the  Dirac  delta-function,  and  the  parameters  P  and 
0  define  the  position  and  the  orientation  of  the  ray.  The 
interferograms  provide  data  on  the  optical  pathlength, 
i( P,0,zq)  ,  along  the  ray  S.  The  reconstruction  problem  is 
to  obtain  the  map  of  n(x,y,zQ)  from  a  number  of  angular  pro¬ 
jections.  A  common  procedure  for  obtaining  the  reconstructed 
map  is  matrix  inversion.  A  square  grid  is  superimposed  over 
the  density  field  and  the  refractive  index  within  each  ele¬ 
ment  (pixel)  is  asstned  to  be  constant.  The  Radon  transfor¬ 
mation  reduces  to  a  set  of  discrete  equations: 


m  m 

)  y  w  (P.0)  n  -  l(P,0)  (2) 

1-1  J-l  J  J 

The  weight  factor  .,  w.  .  are  determined  from  geometric  rela- 

1  J  * 

tions.  Equation  (2)  is  evaluated  for  N  independent  rays, 
resulting  in  N  linear  equations.  N  should  be  larger  or  at 
least  equal  to  the  nmnber  of  pixels,  P.  The  overdeterminate 
system  of  linear  algebraic  equations  can  be  generated  and,  in 


principle,  solved  for  the  unknown  coefficients,  nij(z0)* 


Tomographic  reconstruction  of  the  phase  object  was  achieved 


with  a  modified  version  of  a  well  documented  Algebraic  Recon¬ 
struction  Technique  (ART).  The  linear  equations  (equation  2) 
were  iteratively  solved.  Method  of  conjugate  gradient^  was 
used  to  maximize  the  rate  of  convergance.  A  more  cetailed 
account  of  the  technique  is  given  in  Reference  5. 

Data  Reduction 


A  simplifying  and  common  assunption  made  throughout  this  pro¬ 
gram  has  been  that  the  effects  of  refractive  index  gradient 
on  the  beam  propagation  is  negligible.  This  allowed  the  use 
of  two  dimensional  reconstruction  codes  for  solving  for  the 
phase  map  along  horizontal  planes.  The  reconstructed  data 
consisted  of  the  density  field  at  any  one  elevation  above  the 
airfoil  chord  line.  For  perfect  gas,  and  assuming  constant 
static  pressures  for  each  side  of  the  shockwave,  the  density 
data  may  be  reduced  to  a  velocity  or  mach  number  map. 

Due  to  the  courseness  of  the  calculation  grids  (20  x  20,  or 
100  x  ! 00) ,  the  position  and  the  strength  of  the  shock  may 
not  directly  be  deduced  from  the  tomographic  reconstruction 
results.  A  method  was  devel  ;>ed  to  extract  shock  position 
and  shock  strength  directly  from  the  discontinuities  observed 
in  the  interferograms. 


-5- 


87-2284-13/64 


The  reconstructed  density  field  around  the  leading  edge  of 
the  rotating  blade  at  different  heights  above  the  rotor 
blade,  however,  showed  the  existence  of  a  shock  wave 
extending  beyond  the  blade^.  Determination  of  the  location 
and  the  strength  of  the  shock  wave  using  the  existing  codes 
were  not  possible.  Radon  transformation  in  theory  can  only 
be  applied  to  a  continuous  field  .  Using  an  adaptive 
selection  technique  on  the  original  interferogram  data,  it 
was  possible  to  reconstruct  the  discontinuities  in  the  phase 
map. 

It  was  noted  that  a  number  of  the  interferograms  showed  dis¬ 
continuous  fringe  patterns,  representing  the  passage  of  rays 
through  the  shock  wave.  The  regions  of  the  interferogram 
where  the  fringe  number  distribution  underwent  a  rapid  change 
constituted  the  input  data  for  the  reconstruction  of  the 
shock  waves. 

Figure  4  shows  the  fringe  number  distribution  across  a  shock 
wave.  The  data  of  Figure  2  were  obtained  by  generating  a 
high  resolution  Image  (150  pixels  per  inch)  of  the  fringe 
discontinuity  region.  In  contract,  the  reconstruction  of  the 
entire  field  of  interest  were  achieved  with  a  resolution  of 
10  pixels  per  inch. 
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Reconstruction  of  the  shape  of  the  shock  wave  and  the  shock 


strength  was  based  on  the  assumption  of  Isotropy  for  the  two 
regi  ms  before  and  after  the  shock  wave  (denoted  as  the  ' 
and  regions).  TTie  relation  between  the  flow  Mach  number, 

M,  and  the  Index  of  refraction,  N  -  n  -  can  be  written  as 

M  / M  -  (  1  -  0  N  ' M  2  )  '  '  1  -  (  r-  1  )  C  ,  N'  1  1  '  2  (  i) 

e  ■»  l  m  1  ♦ 

where  Cj  Is  a  constant.  "Hae  normal  shock  assumptions  and 
Giadstone-Da le  formula  relating  the  refractive  index  with  the 
gas  densitv  ran  be  used  to  relate  the  upstream  Mach  number  to 
the  refractive  indices  across  the  shock  wave. 


i  -  c, s  by  2 

l  -  ♦ 

f  4 

1  -  C  N,  M/  ♦  5 

Discontinuous  fringe  patterns  may 

only 

be  ach  ieved  :  or 

ravs 

parallel 

to  the  shock  wave  f  > r  any 

f  i  n  l 

te  lengt  h.  Th  e 

opt  1  - 

cal  p  a  t  h 

difference  between  ravs  1 

a  nd 

2  is  determine) 

bv 

C  (  N 

-  N  )  =  _\F 

(  ■) , 

2  S 

- 

where  I.„ 
s 

is  the  length  )  r  the  she 

k  W  1  '• 

’  >■  el  nme  n  t  pari. 

•  l  e 

the  r  ivs 

,  *.f  Is  the  fringe  o'riber 

1  i  t  t  e 

■  r e n c  e  1 1  ’  r  :  "  t 

!io- 

r on t i nu i 

tv  regions  (see  Figure  1 )  , 

a  nd 

G  ,  is  a  >'  > n s t  a • 

t  . 

An  Iterative  algorithm  was  developed.  Four  unknowns,  M+> 

N_ ,  and  Lg  for  each  "shock  wave  zone”  was  calculated  to 
satisfy  the  equations  (3)  to  (5)  and  a  fourth  condition  on 
the  continuity  of  refractive  indices  around  the  shock  wave  in 
both  regions. 
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2.0  SUMMARY  OF  THE  MOST  IMPORTANT  RESULTS 


a)  Numerical  Experimentation: 

To  evaluate  the  developed  tomographic  reconstruction  code,  a 
numerically  generated  density  ratio  field  shown  in  Figure  5.0 
was  used  as  an  input  for  the  numerical  test.  It  represents 
the  calculated  density  distribution  over  the  tip  region  of 
the  rotor  blade  in  a  plane  above  the  blade^.  The  data  of 
integral  optical  length  difference  as  a  function  of  ray  par¬ 
ameter  9  and  P  was  generated  from  Equation  (2).  The  errors 
associated  with  the  optical  length  measurement,  encountered 
In  the  real  tomography,  were  artificially  added  to  the  data 
file.  Fringe  round-off  errors  of  ±0.5,  ±0.05  and  ±0.005  were 
used  here.  Parametric  study  of  the  reconstruction  code  was 
achieved  for  different  ranges  of  view  angle  and  equation  num¬ 
bers  . 

The  numerical  test  matrix  is  given  in  Table  1.  Figure  6 

shows  the  development  of  the  reconstructed  field  xjj  >rec  w*th 
iteration  number  for  Case  1.  e^,  ig  defined  as  the  average  of 

the  square  of  the  inversion  errors. 

Convergence  rate  of  the  algorithm  and  the  corresponding 
reconstruction  errors  as  functions  of  input  data  errors,  the 
projection  ray  numbers  and  the  view  angles  were  investigated 

-9-  87-2284-1 j/64 
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and  the  results  are  given  in  Figure  7.  The  following  con¬ 


clusions  may  be  drawn  from  Che  numerical  experimentation: 

The  convergence  rate  is  increased  as  the  projection  ray 
nunber,  N,  is  increased  (Figure  7a).  The  corresponding 
inversion  error  for  sufficient  nunber  of  iterations  is,  in 
general,  independent  of  N.  (Figure  7d). 

The  norm  of  the  residual  vector,  R,  in  all  cases,  uncon¬ 
ditionally  converges  to  a  limit  which  is  a  strong  function 
of  the  input  data  error  (Figure  7b).  The  resulting  error 
for  the  first  few  iterations  is  independent  of  the  input 
data  error  (Figure  7e).  However,  for  more  accurate  input 
data,  additional  Iterations  result  in  more  accurate  solu¬ 
tions.  This  underlines  the  importance  of  the  accuracy  of 
the  data  set.  It  also  proposes  a  guideline  for  deter¬ 
mining  the  optimum  nunber  of  iterations  based  on  the  rate 
of  change  of  the  residual. 


Finally,  limited  viewing  angles  (>  3u° )  do  not  have  a  sig¬ 
nificant  effect  on  the  convergence  rate  or  the  accuracy  of 
the  reconstructed  data  (Figures  7e  and  7f).  Acceptable 


angles  were  taken  to  be  symmetric  about  the  axis  of  the 
blade  (90°  to  the  chord  line). 

b)  Reconstruction  of  the  density  and  velocity  field. 

The  reconstructed  density  field  around  the  leading  edge  of 
the  rotating  blade  for  heights  z  »  1.3,  2.5,  3.8  and  5.1  cm 
for  a  7.5  cm  chord  and  blade  aspect  ratio  of  13.7  is  shown  in 
Figure  8.  In  each  case  a  20  x  20  pixel  with  dimensions  of 
10  mm  by  20  mm  was  used.  The  programs  were  run  on  an  IBM  PC. 
The  run  time  for  each  iteration  was  approximately  3  minutes, 
and  the  results  are  shown  after  3  to  6  iterations. 

The  results  show  the  existence  of  the  shock  wave  extending 
beyond  the  blade.  The  magnitude  of  the  maximum  density  ratio 
represents  the  strength  of  the  shock  wave  for  increasing 
height  from  the  blade  surface.  More  detailed  evaluation  of 
the  flow  field  requires  more  refined  mesh  dimensions.  The 
reconstructed  data  at  z  ■  1.3  and  7.8  cm  exhibit  less  noise 
than  the  middle  stations.  This  may  be  attributed  to  the 
error  in  the  input  data. 

The  distribution  of  the  density  ratio  at  three  longitudinal 
positions  x/B  =  0.94,  0.98,  and  1.02  are  shown  in  Figure  9. 
The  data  in  all  cases  exhibit  less  than  adequate  sharpness  in 
the  position  of  the  shock  wave.  This  is  also  attributed  to 
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the  courseness  of  the  reconstruction  grid.  Due  to  limited 
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memory  map,  further  refinement  of  the  grid  on  the  same  micro¬ 
processor  was  not  possible. 

Figure  10  shows  similar  results  in  a  form  of  constant  density 
lines  for  z  -  0.5".  To  transform  the  density  field  into  the 
velocity  field,  accurate  location  of  the  shock  wave  was 
needed  which  was  determined  from  the  fringe  discontinuity 
locations.  Figure  11  shows  the  map  of  constant  velocity  con¬ 
tours.  The  location  of  the  shock  wave  is  shown  as  the  line 
of  discontinuity  in  the  velocity  contours. 
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TABLE  I.  THE  NUMERICAL  TEST  MATRIX 


FRINGE  NUMBER 
ROUND  OFF  ERROR 

VIEW  ANGLE 
(degree) 

NUMBER  OF 
VIEWS 

NUMBER  OF 
POINTS  PER  VIEW 

0.005 

0  ..  180 

61 

20 

0.005 

0  ..  90 

61 

20 

0.005 

0  ..  180 

31 

20 

0.050 

0  ..  180 

61 

20 

0.050 

0  ..  180 

31 

20 

0.500 

0  ..  180 

61 

20 

0.500 

0  ..  180 

31 

20 

0.005 

-75  ..  75 

31 

20 

0.005 

-60  ..  60 

31 

20 

0.005 

-45  ..  45 

31 

20 

0.005 

-30  ..  30 

31 

20 

0.005 

-15  ..  15 

31 

20 

0.050 

-75  ..  75 

31 

20 

0.050 

-60  ..  60 

31 

20 

0.050 

-45  ..  45 

31 

20 

0.050 

-30  ..  30 

31 

20 

0.050 

-15  ..  15 

31 

20 
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Figure  5.  Numercially  generated  density  and  refractive 
index  distribution  (Ref,  7), 
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Figure  II.  Reconstructed  Velocity  Contours 


